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ABSTRACT 

MR image sparsity/compressibility has been widely exploited 
for imaging acceleration with the development of compressed 
sensing. A sparsity-based approach to rigid-body motion cor- 
rection is presented for the first time in this paper. A motion 
is sought after such that the compensated MR image is max- 
imally sparse/compressible among the infinite candidates. It- 
erative algorithms are proposed that jointly estimate the mo- 
tion and the image content. The proposed method has a lot 
of merits, such as no need of additional data and loose re- 
quirement for the sampling sequence. Promising results are 
presented to demonstrate its performance. 

1. INTRODUCTION 

Imaging artifacts due to patient motion remain a major chal- 
lenge in many MRI applications. In this paper, we consider 
MR image reconstruction from fc-space data corrupted by 2D 
translational motion. The problem is ill-posed since there 
exist infinitely many image candidates and possible motions 
which result in the same fc-space data. The rigid-body mo- 
tion correction problem has been an active research topic for 
a long time. Existing methods resolve the problem by acquir- 
ing more fc-space data and/or adopting a specific sampling 
sequence. Early works assume that the region of in- 

terest (ROI) is a priori known and do motion correction by 
borrowing the idea of phase retrieval [3 |. A relatively small 
ROI is helpful for the motion correction and is obtained by 
oversampling in the frequency domain, or equivalently, im- 
posing an enlarged field of view (FOV) in the image domain. 
Navigator-based methods |4,5 | have been popular in the past 
decade which resolve the ill-posed problem by acquiring ex- 
tra fc-space lines named as "navigator". The motion detected 
between the navigator and a motion-free reference is used 
for data correction. Translational motion is studied in Q 
while |5| also considers rotational motion. In 16), the mo- 
tion is estimated by exploiting correlations between adjacent 
readout lines by assuming that little motion exists during a 
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single echotrain (the same assumption required in navigator- 
based methods). To guarantee the existence of the correla- 
tions, oversampling is required along the phase-encode direc- 
tion ||6). 

Compressed sensing (CS) fl] has recently become a pow- 
erful tool for signal processing that aims at reconstructing 
a high dimensional signal from its low dimensional linear 
measurements. The inherent reason of its success is that 
most natural signals can be sparsely represented in some ba- 
sis, rendering that a small number of Unear measurements 
are sufficient for the signal recovery. Since medical images 
are sparse/compressible, for example, in a wavelet domain, 
applications of CS to MRI have been vastly reported in the 
literature (e.g., (W\) for accelerating the scanning process by 
downsampling in the fc-space. 

In this paper, we study the motion correction problem by 
exploiting the MR image sparsity. To do this, we formulate 
the motions as unknown parameters in the sampling system. 
Rather than acquiring partial /c-space data in the motion-free 
environment for imaging acceleration, we attempt to estimate 
the system parameters jointly with the image reconstruction 
process by exploiting the redundancy in the full-sampled fc- 
space data. The idea is inspired by our recent work |9 1 where 
it is shown that uncertain system parameters can be accurately 
estimated along with the sparse signal in one particular situ- 
ation. The proposed method relies on the assumption that 
little motion occurs during a single readout line, which can 
be approximately satisfied in practice and is greatly relaxed 
in comparison with that in [4,6]. Unlike existing navigator- 
based methods, it does not require additional fc-space data. 
The method is justified using a Shepp-Logan phantom and 
simulated human brain data. 

2. PROBLEM FORMULATION OF MOTION 
CORRECTION 

Consider that a translational motion f3^ occurs when acquir- 
ing the fc-space measurement at k. So it holds at the moment 
m (r) — TO° (r — /3^), where m° is the MR image of inter- 
est, m is the translated image and r denotes the coordinate in 



the image domain. According to the relationship between an 
MR image and its fc-space measurements, it holds 
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where M° — Fm° and M denote the motion-free and 
motion-corrupted fc-space data, respectively, with F denot- 
ing the 2D Fourier transform operator. It is shown in ([TJ 
that the translational motion leads to a phase error and un- 
altered amplitude. To represent ([T]i more compactly, let 



/3 



and then write 



mto 
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where Tp denotes a linear operator caused by the translational 
motion (3 and is referred to as the translational operator here- 
after. In particular, let Ap be a matrix of the same dimen- 
sion as M with its element kp (fc) = e"*^''^'"' '^'"^ Then 
7^M° = M° where denotes the Hadamard prod- 



uct. 



In addition, it holds that Tg 



= T-/3 and To =1 where 
I denotes the identity operator. For a particular motion /3, the 
linear system of equations in Q relates the MR image mf of 
interest and the acquired fc-space data M. Our objective is 
to reconstruct the image m° and possibly the motion j3 given 
the corrupted fc-space data M. 

The challenge is that the problem of recovering m° from 
M is ill-posed since there exist infinite number of solutions 
to (pji. In particular, for any (3 there exists a solution m = 
F T-pM such that (j2| holds. To choose the correct one 
among the infinite candidates, as a result, additional informa- 
tion has to be exploited. 

3. SPARSITY-DRIVEN MOTION CORRECTION 
3.1. Sparsity-based formulation 

We first consider that the MR image is directly constructed 
from the corrupted fc-space data M without considering the 
translational motion. Then, the constructed image is 

m = F-^M = T-^ (A^ M°) = (j-^A^g) 0m°, (3) 

where denotes the circular convolution operation. So the 
obtained image is the true image ra° after a convolution with 
Ap. That explains how the imaging artifacts come from 
the translational motion. Due to the imaging artifacts, it is nat- 
ural to conjecture that the translational motion will reduce the 
sparsity/compressibility of MR images, i.e., the true image is 
the maximally sparse/compressible solution (under an appro- 
priate basis) to Q. Examples of a Shepp-Logan phantom and 
simulated human brain are presented in Fig. [T] In compari- 
son with the motion-free images (in col 1), severe artifacts are 
present in the motion-corrupted ones (in col 2). Numerically, 



it can be shown that the l\ norms (a commonly used spar- 
sity metric) of the motion-corrupted images (under an Haar 
wavelet basis) are larger than those of the motion-free ones. 

Based on the conjecture above, we propose to reconstruct 
the motion-free MR image using the maximally sparse solu- 
tion that satisfies the data consistency constraint (j2]i, i.e., by 
solving a basis pursuit (BP) like optimization problem: 



subject to TpTm = M and f3 G Dp, 



mm 
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where >V is a sparsifying operator (e.g., a wavelet transform), 
||-||-^ denotes the £i norm (sum of amplitude of all elements) 
that is commonly used to promote sparsity and Dp denotes 
the domain of f3. Given an appropriate constant C, an equiv- 
alent (LASSO-like) formulation of (HI is 



min \\M — TpJFmW , 

m,j3 ^ 

subject to ||>Vm||^ < C and /3 e Vp, 



(5) 



where 1 1 • 1 1 p denotes the Frobenius norm and C" < 1 1 ^ M 1 1 
is a constant that needs to be tuned in practice and can be set 
to C = ||>Vm°||j in an ideal case. The motions P^. with 
respect to the coordinate k are not independent. For example, 
the motion should be piecewise smooth when sorted chrono- 
logically. In this paper, we assume that the same motion 
is shared among every readout line, which can be approxi- 
mately satisfied in practice and greatly reduces the number 
of unknown parameters. Note that the assumption is greatly 
relaxed comparing to that in |4,6| where it is assumed that 
little motion exists during every echotrain (each compris- 
ing several readout lines). Moreover, /3 should be properly 
bounded in practice. These a priori knowledge narrows the 
selection of f3 and defines Dp. To the best of our knowl- 
edge, Q and (|5]) present the first sparsity-based formulations 
for the MRI motion correction problem though CS has been 
widely applied to MRI for imaging acceleration. 

3.2. An alternating algorithm 

Due to the nonconvexity with respect to f3, problem (|5]l is 
nonconvex. We propose an alternating algorithm that itera- 
tively carries out the following three steps starting with j = 
and /3(") = 0: 



1) solving 



m 



0+1) 



:argmin llM - Tg(j)Tm\\ . 

m ^ "F 

subject to ||>Vm||^ < C. 



(6) 



2) solving 



M - TpTm 



(j+i) 



(7) 



subject to /3 e Vp. 



3) letting J ^j + l. 

To interpret the algorithm above, denote two sets 5*1 = 
{m : \\Wm\\^ < CjandSa = {m :M = TpTm,f3 e Vp}. 
We refer to 5*1 as the sparse domain and 5*2 as the Fourier 
domain, respectively. Denote Vi and V2 projections onto 5*1 
and S2, respectively. Let m'-'' = jF~^T_0U)M ■ Then it 



holds m 



Vim^^^ following from Step 1. At Step 2, 



the equality m(^+^^ = 7^2 'ti^"'^^'' holds. So the algorithm 
above is equivalent to the recursion 



V2Vim 



{3} 



(8) 



starting with m^"^ — rh where m is as defined in (jsjl. Since 
m'") e 5*2, the recursion attempts to find a point in ^2 near- 
est to Si, or equivalently, the maximally sparse image that is 
consistent with the fc-space observation. The constant C is 
generally unavailable in advance and needs to be tuned in the 
algorithm. A simple method is to set C such that the obtained 
image has the least ^1 norm in the sparsifying domain (which 
is consistent with our objective). Vi (i.e., to solve (|6]l) is easy 
to compute. The next subsection is devoted to computing 'P2- 

3.3. Approximate computation of V2 

The computation of V2 is equivalent to solving the optimiza- 
tion problem in (j?]). Under the assumption that a common 
displacement is shared among each readout line and the mo- 
tions among different readout lines are independent, the num- 
ber of unknown variables in /3 is reduced to twice the number 
of readout lines (each displacement has 2 degrees of freedom 
in the 2D case) and the displacement in each line can be sep- 
arately calculated. However, the exact calculation of V2, or 
equivalently, to obtain the global optimum of (|7]l, is generally 
infeasible since the Fourier domain S2 is severely nonconvex. 
The /3 solution is typically trapped at a local minimum if op- 
timization methods, e.g., a gradient method, are used to solve 
(|7]i. Inspired by a navigator-based method in |5 1, we propose 
a practically efficient algorithm as follows. 

In the term Tm^-'\ denoted by M^^\ can be consid- 
ered as the current estimate of the motion-free fc-space data 
M° (without accounting for the difference between m*^"''' and 
m (j)). To solve ^ is in fact to estimate the translational mo- 
tion between the motion-corrupted observation M and the 
current motion-free estimate M^-'K So we may consider each 
readout line of M as a navigator and the associated line of 
M'^-'^ as the reference. Then the motion within each readout 
line can be estimated using the navigator-based method in |j5). 
We omit the details. At last, we point out that the algorithm 
performance can be improved in practice by a modification 
of the current motion-free fc-space estimate: Af*--' -' = | Af | 

sgn ^J^m^"'''^, instead of using Af*^-'-' — J^rh^-'\ where the 

absolute operator | • | does an element-wise operation. The un- 
derlying reason is obvious according to ([T]). The translational 



motion changes only the phase information of the fc-space 
data and keeps the amplitude. So intuitively, the modified es- 
timate is closer to the true motion-free data and leads to better 
performance. 

3.4. Sparse RAAR for motion correction 



According to ([8]), the algorithm in Subsection 3.1 is imple- 
mented by iteratively projecting onto the two sets 5*1 and 5*2. 
It is related to iterative projection algorithms {3 \ for phase re- 
trieval, which aims to find an intersection point of two sets 
and are built upon combining projections onto the two sets in 
some fashion. In fact, the recursion ([8]l is exactly the error 
reduction (ER) algorithm 1 3J without accounting for the dif- 
ferences of the sets. One drawback of ER is its slow conver- 
gence. In our setting, we want to solve a feasibility problem: 



find m e S"! n 52, 



(9) 



i.e., to find a point lying in both the sparse domain Si and the 
Fourier domain ^2. The relaxed averaged alternating reflec- 
tions (RAAR) algorithm introduced in ^OJ has fast conver- 
gence speed and stable performance. RAAR starts with some 
initial point m^'^^ and is defined by the recursion 



(7^l7^2+I) + (l 



(10) 



where 9 € [0, 1] is a constant, and reflectors TZi = 27^1 — I 
and Tl2 = 27^2 — I- By defining and computing Vi and 7^2 
as before, we propose an algorithm for the motion correction 
problem, named as sparse RAAR (SRAAR), where m in ([3]l 
is used as the initial point. 

4. RESULTS 

The proposed method is validated using a Shepp-Logan phan- 
tom and simulated human brain data obtained from BrainWeb 
[^with both of size 256 x 256. Continuously varying transla- 
tional motions among the readout lines are randomly gener- 
ated. Motion artifacts are introduced by applying varying lin- 
ear phase shifts to the motion-free fc-space data according to 
([T]l. An Haar wavelet is selected as the sparsifying transform. 
SRAAR is applied to the motion-corrupted fc-space data to 
reconstruct the MR images with the setting 6 — 0.9. SRAAR 
is terminated after a fixed number of iterations. 

Simulation results are presented in Fig. [T] For the phan- 
tom (row 1), almost exact reconstruction is obtained using 
SRAAR. A small amount of motions (within 5 pixels along 
both the readout and phase-encode directions) are studied in 
rows 2 (without noise) and 3 (noise added) with the human 
brain. It can be seen that even a small amount of motions 
may cause severe imaging artifacts. After the sparsity-based 
motion correction with SRAAR, only few artifacts remain. In 



'http://www.bic.mni.mcgill.ca/brainweb 



Motion-free Motion-corrupted SRAAR reconstruction Difference 




Fig. 1. Simulation results of sparsity-based motion correc- 
tion on a Shepp-Logan phantom (row 1, 100 iterations), sim- 
ulated human brain (row 2, 200 iterations), and noisy brain 
with small motions (row 3, 200 iterations) and large motions 
(row 4, 400 iterations). 

the presence of both noises and a large amount of motions (as 
3 times large as those in rows 2 and 3), it is shown in row 4 
that SRAAR may have difficulties to produce a good result 
though most artifacts are removed. We note that SRAAR is 
computationally efficient in general. In the above, each iter- 
ation takes about Is in Matlab v.7.7.0 on a PC with a 3GHz 
CPU, i.e., each reconstructed image is obtained within few 
minutes. Moreover, SRAAR can be greatly accelerated by 
estimating motions in the readout lines in parallel when com- 
puting 7^2. 

5. CONCLUSION 

A first sparsity-based approach to motion correction was in- 
troduced in this paper. An efficient algorithm was proposed 
partially inspired by phase retrieval and existing navigator- 
based methods. The promising results presented in the paper 
indicate good potential for practical application of the pro- 
posed method. The current work considers only translational 
motion and encourages further studies of more complicated 
motions. Reference |9| shows that unknown system parame- 
ters can be estimated even in the undersampling case, suggest- 
ing that imaging acceleration is also possible in the presence 
of motions. A shortcoming of the algorithms in this paper is 



the need of the tuning of C. A recent improvement in the 
algorithm design can be found in pT| . 
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